#This script executes event-study regressions to estimate cohort-specific 
#associations between first-generation prenatal exposure to nonattainment 
#and first-generation family structure outcomes.

rm(list=ls())
gc()
library(reldist)
library(ineq)
library(stargazer)
library(data.table)
library(dplyr)
library(dtplyr)
library(ggplot2)
library(stringdist)
library(lfe)
library(haven)
library(readr)
library(broom)
library(tidylog)
library(rdrobust)

# Set the working directory
setwd("/projects/opp_env/caa1970/")

# Source external R scripts for custom functions
source("Code/Child_Outcomes/child_outcomes_functions.R")
source("Code/rounding.r")

# Load the first-generation analysis dataset with ACS link.
load("Data/first_gen_analysis_data_jepmicro_acs.rda")

# Load the first-generation analysis dataset
load("Data/first_gen_analysis_data_jepmicro.rda")

# Event Studies

#Any Links

eventB3a <- felm(any_link ~inst_1969+inst_1970+inst_1972+
                 inst_1973+inst_1974+inst_1975 |state_year+county_factor+year_factor+month+
                 survey_year_by_year | 0 |county_factor,
                 data=parent_piks%>% filter(year%in%1969:1975, AGE >33 & AGE < 44, real_wages < wage_cap, !is.na(all_fullterm), !is.na(lfpr), !is.na(PI_PC)))
summary(eventB3a)

#Save the underlying point
figB3a <- eventB3a %>% tidy() %>% mutate_each(funs(signif(.,4)),-term)  %>% 
  filter(grepl("inst", term), !is.na(estimate))%>%
  mutate(year=as.numeric(gsub("inst_","", gsub("_bin", "", term)))) %>% 
  bind_rows(., data.frame(estimate=0, std.error=0, year=1971)) %>% filter(year<=1975) 
write.csv(figB3a <- eventB3a %>% tidy() %>% mutate_each(funs(signif(.,4)),-term)  %>% 
, file="JPE_Micro/Output/B3_Figures/FigB3a.csv", row.names=F)

# Number of Links

eventB3b <- felm(total_any_links~inst_1969+inst_1970+inst_1972+
                 inst_1973+inst_1974+inst_1975 |state_year+county_factor+year_factor+month+
                 survey_year_by_year | 0 |county_factor,                
               data=parent_piks%>% filter(year%in%1969:1975, AGE >33 & AGE < 44, real_wages < wage_cap, !is.na(all_fullterm), !is.na(lfpr), !is.na(PI_PC)))
summary(eventB3b)

#Save the underlying point
figB3b <- eventB3b %>% tidy() %>% mutate_each(funs(signif(.,4)),-term)  %>% 
  filter(grepl("inst", term), !is.na(estimate))%>%
  mutate(year=as.numeric(gsub("inst_","", gsub("_bin", "", term)))) %>% 
  bind_rows(., data.frame(estimate=0, std.error=0, year=1971)) %>% filter(year<=1975) 
write.csv(figB3b, file="JPE_Micro/Output/B3_Figures/FigB3b.csv", row.names=F)

## Second Gen Sample Link

eventB3c <- felm(second_gen~inst_1969+inst_1970+inst_1972+
                 inst_1973+inst_1974+inst_1975 |state_year+county_factor+year_factor+month+
                 survey_year_by_year | 0 |county_factor, 
               data=parent_piks_acs%>% filter(year%in%1969:1975,AGE >33 & AGE < 44, real_wages < wage_cap, !is.na(all_fullterm),!is.na(parent_age_min), !is.na(lfpr), !is.na(PI_PC)))
summary(eventB3c)

#Save the underlying point
figB3c <- eventB3c %>% tidy() %>% mutate_each(funs(signif(.,4)),-term)  %>% 
  filter(grepl("inst", term), !is.na(estimate))%>%
  mutate(year=as.numeric(gsub("inst_","", gsub("_bin", "", term)))) %>% 
  bind_rows(., data.frame(estimate=0, std.error=0, year=1971)) %>% filter(year<=1975) 
write.csv(figB3c, file="JPE_Micro/Output/B3_Figures/FigB3c.csv", row.names=F)

## Age at First Link

eventB3d <- felm(parent_age_min~inst_1969+inst_1970+inst_1972+
                 inst_1973+inst_1974+inst_1975 |state_year+county_factor+year_factor+month+
                 survey_year_by_year | 0 |county_factor,                 
               data=parent_piks_acs%>% filter(year%in%1969:1975, AGE >33 & AGE < 44, real_wages < wage_cap, !is.na(all_fullterm), !is.na(lfpr), !is.na(PI_PC)))
summary(eventB3d)

#Save the underlying point
figB3d <- eventB3d %>% tidy() %>% mutate_each(funs(signif(.,4)),-term)  %>% 
  filter(grepl("inst", term), !is.na(estimate))%>%
  mutate(year=as.numeric(gsub("inst_","", gsub("_bin", "", term)))) %>% 
  bind_rows(., data.frame(estimate=0, std.error=0, year=1971)) %>% filter(year<=1975) 
write.csv(figB3d, file="JPE_Micro/Output/B3_Figures/FigB3d.csv", row.names=F)

## Married with kids

modelB3e <- felm(married~inst_1969+inst_1970+inst_1972+
                 inst_1973+inst_1974+inst_1975 |state_year+county_factor+year_factor+month+
                 survey_year_by_year | 0 |county_factor,               
               data=parent_piks_acs%>% filter(year%in%1969:1975, AGE >33 & AGE < 44, real_wages < wage_cap, !is.na(all_fullterm), !is.na(lfpr), !is.na(PI_PC)))
summary(modelB3e)

#Save the underlying point
figB3e <- modelB3e %>% tidy() %>% mutate_each(funs(signif(.,4)),-term)  %>% 
  filter(grepl("inst", term), !is.na(estimate))%>%
  mutate(year=as.numeric(gsub("inst_","", gsub("_bin", "", term)))) %>% 
  bind_rows(., data.frame(estimate=0, std.error=0, year=1971)) %>% filter(year<=1975) 
write.csv(figB3e, file="JPE_Micro/Output/B3_Figures/FigB3e.csv", row.names=F)

## Divorced with kids

modelB3f <- felm(divorced~inst_1969+inst_1970+inst_1972+
                 inst_1973+inst_1974+inst_1975 |state_year+county_factor+year_factor+month+
                 survey_year_by_year | 0 |county_factor,                
               data=parent_piks_acs%>% filter(year%in%1969:1975, AGE >33 & AGE < 44, real_wages < wage_cap, !is.na(all_fullterm), !is.na(lfpr), !is.na(PI_PC)))
summary(modelB3f)

#Save the underlying point
figB3f <- modelB3f %>% tidy() %>% mutate_each(funs(signif(.,4)),-term)  %>% 
  filter(grepl("inst", term), !is.na(estimate))%>%
  mutate(year=as.numeric(gsub("inst_","", gsub("_bin", "", term)))) %>% 
  bind_rows(., data.frame(estimate=0, std.error=0, year=1971)) %>% filter(year<=1975) 
write.csv(figB3f, file="JPE_Micro/Output/B3_Figures/FigB3f.csv", row.names=F)